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ABSTRACT 


A  study  is  made  of  the  effect  of  a  nonlinearity  on  the  stability 
properties  of  a  simplified  version  of  the  PIC  difference  equations*  An 
equilibrium  amplitude  of  the  fluctuations,  resulting  from  the  instability, 
is  determined  as  a  function  of  the  parameters  of  the  method. 
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INTRODUCTION 


The  Particle-in-Cell  (PIC)  technique  is  a  finite  difference  method 
of  expressing  the  equations  of  motion  of  a  compressible  fluid;  it  has 
been  described  in  various  reports  [1,2].  The  computational  framework 
for  the  method  is  achieved  by  dividing  the  system  into  an  Eulerian  mesh 
of  cells  and  superimposing  a  mesh  of  particles  whose  distribution  and 
mass  are  such  as  to  describe  the  initial  configuration  of  the  fluid. 

The  differential  equations  of  motion,  with  transport  terms  neglected, 
are  written  in  finite  difference  form  relative  to  the  system  of  cells. 

The  transport  effect  is  obtained  by  allowing  the  particles  to  move 
through  the  fixed  cellular  coordinate  system  according  to  the  velocities 
in  the  neighboring  cells. 

The  method  has  achieved  considerable  success  in  problems  involving 
high  velocity  fluid  flow  and  it  has  been  recognized  [ 1 J  that  a  large 
portion  of  this  success  is  attributable  to  the  transport  mechanism  of 
the  method.  The  movement  of  particles  across  cell  boundaries  gives  rise 
to  a  nonlinear  dissipative  force  which  is  effective  in  reducing  the  fluc¬ 
tuations  that  arise  as  a  result  of  the  differencing  technique.  This  dis¬ 
sipative  term  is  of  the  form  of  a  "true"  viscosity,  in  that  it  is  propor¬ 
tional  to  the  velocity  gradient. 
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The  same  measure  of  success  has  not  been  obtained,  however,  in  rep¬ 
resenting  low  velocity  flow.  The  primary  reason  for  this  is  that  the 
velocity  gradients  are  too  small  to  make  the  dissipative  force  effective, 
so  that  instabilities  develop.  This  difficulty  can  be  overcome  by  intro¬ 
ducing  into  the  equations  linear  forms  of  artificial  viscosity,  whose 
effect  persists  even  to  zero  speeds.  However,  this  solution  is  not 
always  desirable  in  other  regions  of  the  fluid,  since  the  artificial  vis¬ 
cosity  tends  to  obliterate  some  of  the  features  of  interest  of  the  high 
velocity  flow. 

But  the  effect  of  this  instability  can  be  minimized,  without  re¬ 
course  to  artificial  viscosity,  by  the  optimum  choice  of  the  adjustable 
parameters  of  the  method.  The  reason  is  that  the  dissipative  term  pre¬ 
vents  unbounded  growth  of  the  instability.  Hence,  if  one  can  determine 
the  upper  limit  of  the  fluctuations,  which  result  from  this  instability, 
as  a  function  of  the  parameters  of  the  system,  the  parameters  can  then 
be  chosen  in  such  a  way  as  to  bound  the  instability  to  a  tolerable  level. 
The  primary  purpose  of  this  report  is  to  determine  the  mechanics  of  this 
dissipative  process  in  a  finite  difference  system  and  thereby  achieve  an 
expression  for  the  equilibrium  amplitude  of  fluctuations  in  terms  of  the 
parameters. 

This  study  was  further  prompted  by  a  desire  to  learn  more  about  the 
effect  of  nonlinearities  on  difference  equation  stability  in  general. 

This  consideration,  together  with  the  difficulty  involved  in  the  study 
of  more  than  one  nonlinearity  at  a  time,  has  led  us  to  investigate  a 
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simplified  version  of  the  PIC  equations  rather  than  the  full  equations. 

It  is  hoped  that  the  results  of  this  analysis  will  thereby  be  applicable 
to  a  wider  range  of  difference  methods. 

The  appendices  to  this  report  deal  with  certain  improvements  on  the 
PIC  method  which  have  not  previously  been  published.  Appendix  I  is  con¬ 
cerned  with  the  form  of  the  energy  equation  as  it  appears  in  this  report. 
This  form  differs  from  that  which  was  employed  in  previous  discussions  of 
the  PIC  method  in  that  it  is  derived  from  the  differential  equation  of 
motion  for  specific  internal  energy  rather  than  that  for  total  energy. 

Appendix  II  describes  the  form  and  effectiveness  of  the  various  types 
of  artificial  viscosity  which  have  been  employed  in  PIC  calculations  at 
Los  Alamos.  It  also  contains  a  discussion  of  a  more  effective  manner  of 
expressing  the  viscosity  terms  in  the  difference  equations  than  had  been 
in  use  previously.  This  latter  section  should  also  be  applicable  to 
finite  difference  methods  other  than  PIC. 
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PART  I.  STABILITY  ANALYSIS  OF  THE  PIC  EQUATIONS 


A  PIC  calculation  is  performed  in  three  phases.  In  Phase  I  the  dif¬ 
ference  form  of  the  momentum  and  energy  equations,  exclusive  of  the  trans¬ 
port  terms,  is  solved.  Phase  II  is  then  concerned  with  the  movement  of 
particles;  Phase  III  deals  with  the  re-partitioning,  among  the  cells,  of 
momentum  and  energy  carried  by  particles  which  cross  cell  boundaries. 

The  stability  properties  of  the  method  are  completely  determined, 
however,  by  the  first  of  these  phases,  since  they  are  concerned  with  the 
response  of  the  system  to  a  perturbation  of  steady  state  conditions. 

For,  if  the  fluctuations  which  result  from  this  perturbation  attain  suf¬ 
ficient  magnitude  to  cause  appreciable  movement  of  particles,  then  in¬ 
stability  is  amply  demonstrated.  But,  by  restricting  the  study  to  the 
equations  of  Phase  I,  the  determination  of  stability  becomes  a  matter  of 
ascertaining  whether  or  not  there  is  growth  with  time  of  the  magnitude 
of  the  fluctuations  produced  by  the  perturbation. 

Assuming,  therefore,  that  the  system  will  retain  its  initial  uni¬ 
form  density  configuration,  the  one  dimensional  PIC  difference  equations 
can  be  written. 
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where  a  and  f  are  constants  of  approximately  unit  magnitude  and  cQ  is  a 
fixed  representative  sound  speed  for  the  system. 

The  q  terms  represent  artificial  viscosity;  their  purpose  is  to  im¬ 
prove  the  stability  and  accuracy  of  the  method  by  smoothing  out  any  fluc¬ 
tuations  which  may  develop  in  the  system.  A  discussion  of  the  form  of 
these  terms  and  the  manner  in  which  they  are  expressed  in  the  difference 
equations  will  be  found  in  Appendix  II,  but  it  is  pertinent  to  this  sec¬ 
tion  to  point  out  one  particular  aspect  of  their  expression  which  compli¬ 
cates  the  analysis.  Artificial  viscosity  is  usually  applied  only  in 
those  parts  of  the  system  that  are  undergoing  compression;  in  those  re¬ 
gions  which  are  experiencing  rarefaction,  additional  dissipation  is  not 
required. 
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The  dilemma  which  arises  in  a  stability  analysis  is  that  there  is  no 
clear cut  way  of  expressing,  in  the  equations,  the  fact  that  the  q  terms 
are  nonzero  only  when  their  velocity  gradient  factor  is  positive.  The 
course  which  is  sometimes  followed  in  such  a  situation  is  to  assume  that, 
in  a  perturbed  stagnation,  each  point  in  the  system  experiences  equal 
periods  of  compression  and  rarefaction.  Then,  for  the  purpose  of  the 
analysis,  the  dissipative  mechanism  is  applied  at  all  times  but  only  one- 
half  the  actual  coefficient  of  viscosity  is  used.  For  lack  of  an  alter¬ 
native,  this  procedure  will  be  followed  in  the  present  analysis,  but  the 
amount  of  error  thereby  introduced  will  be  demonstrated  by  experiment. 

To  proceed  with  the  stability  analysis,  let  us  consider  a  system, 
described  by  Eqs.  (1),  which  has  been  perturbed  from  steady  state.  The 
perturbation  is  introduced  by  means  of  the  substitution 

u“-i/a  -*uo  +  6Vi/2’ 

where  uQ  and  1^  are  the  steady  state  values  and 


Then,  retaining  only  first  order  terms  in  the  variation  of  these  small 
values,  Eqs.  (1),  with  a  polytropic  equation  of  state,  p  =  (7  -  1)pl,  be. 


come 
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Assuming  that  a  solution  of  Eqs.  (2)  can  be  expressed  in  terms  of 
Fourier  series,  let  us  examine  the  conditions  under  which  a  typical  term 
can  be  a  solution.  Take  as  the  representative  terms  from  the  series 


n  ik(j-l/2)  n 

&u._l/2=  Ae 

RT»  =  B  eik(j-V2)  rn 

5I0-l/2  B  6  r  * 


(3) 


Then  expressing  Eqs.  (2)  in  terms  of  solution  (3),  we  obtain 

,\  -ik  ik,  -ik  ik  , 

A(r  -  1)  =  ffB(e  -  e  )  +  AA(e  +  e  -  2), 


oI0 

B(  r  -  1 )  -  A| 


/  —ik  ik,  t  1  \ 
(e  -  e  )  (r  +  1) 


and,  on  collecting  terms,  this  becomes 

A[r  -  1  +  2A(1  -  cos  k)]  +  B(2icr  sin  k)  =  0, 

w 

A[ioIQ(r  +1)  sin  k]  +  B(r  -  1)  =0. 

A  necessary  and  sufficient  condition  that  Eqs.  (4)  have  a  nontrivial  solu¬ 
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tion  for  A  and  B  is  that 


(r  -  I)2  +  2A(1  -  cos  k)  (r  -  1)  +  2a2IQ  (sin2  k)  (r  +  1)  =  0  . 


From  this  equation  a  solution  can  be  obtained  for  r  in  the  form 


where  a  =  A(l  -  cos  k) , 

2  .  2 

P  =  a  IQ  sin  k. 

In  order  that  solution  (3)  does  not  grow  with  time  it  is  necessary 
that  | r |  <  1.  The  restrictions,  which  this  condition  imposes  upon  the 
parameters,  depend  on  whether  r  is  a  real  or  complex  number. 

2 

Case  I:  (a  +  p)  <  4p 

Then  r  is  complex  and 

|  r  | 2  =  1  -  2a  +  2P 

O  2 

=  1  -  2A(1  -  cos  k)  +  2a  Iq(1  -  cos  k) . 

Thus,  in  order  that  |r|  <  1,  it  is  required  that 

A  >  a2I0( 1  +  cos  k) . 

Now,  the  boundary  conditions  on  a  PIC  problem  require  that  the  ve¬ 
locity  vanish  at  the  ends  of  the  system.  This  limits  the  allowable  values 
of  k  to  the  values 


where  N  is  the  number  of  cells  in  the  system.  Thus,  to  insure  stability 
for  all  values  of  k,  it  is  sufficient  to  require  that  the  inequality 
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above  holds  when  k  =  The  stability  property  for  Case  I  then  becomes 


(7  -  1)  IQ  5t 

*c0  +  f|uj  >  (1  +  cos  -)  - 
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But  at  steady  state  uQ  s  0,  so  for  Case  I  the  stability  requirement  is 


acn  >  0  +  cos  -) 
>2 


K  <r-  i)  i0at 


0  -  '  '  '  N'  2&x 
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Case  II:  (a  +  p)  > 

In  this  case  r  is  real  and  c)r/Sk  is  proportional  to  sin  k,  so  that 

| r |  attains  its  maximum  value  at  k  =  it.  For  this  value  of  k  the  require¬ 
ment  that  | r |  <  1  reduces  to  the  condition 


i.e. 


0  <  2*  <  1 , 


>  ox 

0  <  ao0  &t- 


Notice,  in  Eqs.  (5),  that  k  =  it  corresponds  to  a  solution  with  a  wave¬ 
length  of  2  cells.  Hence  it  is  the  high  frequency  components  which  are 
most  drastically  affected  by  excessive  artificial  viscosity.  The  reason  is 
that  the  viscosity  term  is  proportional  to  the  velocity  gradient  between 
adjacent  cells,  so  that  excessive  correction  results  in  a  gradient  of  oppo¬ 
site  direction  but  increased  magnitude.  The  resultant  of  these  excessive 
corrections  is  a  high  frequency  disturbance  which  grows  exponentially. 


^Fourier  components  for  which  k  =  2nn  correspond  to  velocity  profiles 
which  are  identically  zero;  hence  growth  of  these  components  will  not 
lead  to  instability. 
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The  result  of  violating  the  condition  which  is  derived  in  Case  I, 
i.e.,  too  little  artificial  viscosity,  is  the  subject  of  Part  II. 

Briefly,  rapidly  growing  fluctuations  develop  here  also,  until  the  ve¬ 
locity  magnitudes  become  large  enough  to  allow  the  f |Uq|  term  to  sta¬ 
bilize  the  system.  This  capacity  permits  PIC  calculations  to  be  per¬ 
formed  without  any  artificial  viscosity  whatsoever,  since  a  term  of  the 
f |u0|  type  is  inherent  in  the  PIC  method  of  moving  the  particles  (Ref.  1, 
p.  14  ff). 

A  plot  of  the  predicted  region  of  stability  (for  fixed  values  of  7, 
5x,  N  and  IQ)  is  shown  in  Fig.  1 .  When,  in  the  machine  runs,  artificial 
viscosity  is  applied  in  both  rarefactions  and  compressions,  it  is  found 
that  the  observed  limits  of  stability  agree  perfectly  with  this  predic¬ 
tion.  But  the  experimentally  observed  outline  of  the  region  of  sta¬ 
bility  for  problems  in  which  artificial  viscosity  is  applied  in  compres¬ 
sions  only  is  shown  by  the  dashed  line.  The  discrepancy  between  the  area 
thus  outlined  and  the  predicted  region  of  stability  indicates  that  our 
assumption  (that  in  a  perturbed  stagnation  artificial  viscosity  applied 
in  compressions  only  is  equivalent  to  one-half  the  amount  applied  in  both 
compressions  and  rarefactions)  does  not  hold  for  large  amounts  of  artifi¬ 
cial  viscosity.  The  reason  is  that,  when  compressions  are  rapidly 
damped,  the  system  experiences  rarefaction  considerably  more  than  fifty 
percent  of  the  time. 
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PART  II.  THE  NONLINEAR  BOUNDING  OF  THE  PIC  INSTABILITY 

The  stability  analysis  of  Part  I  indicated  that  the  PIC  equations 
ve re  unstable  when  the  coefficient  of  the  linear  part  of  the  artificial 
viscosity  term  was  less  than  a  certain  positive  definite  function  of  the 
parameters  of  the  system.  It  was  also  observed,  however,  that  this  in¬ 
stability  was  bounded  in  amplitude  by  the  nonlinear  part  of  the  viscosity 
terms,  which  corresponds  to  the  effective  viscosity  of  PIC.  It  is  the  pur¬ 
pose  of  this  section  of  the  report  to  study  the  mechanics  of  this  damping 
force  and  to  obtain  an  expression  for  the  equilibrium  amplitude  of  fluc¬ 
tuations  in  terms  of  the  parameters  of  the  system. 

To  attempt  an  analysis  of  the  complete  PIC  equations,  including  a 
provision  for  density  variation,  would  be  a  very  complicated  procedure. 
Instead,  the  study  will  be  concerned  with  the  effects  of  the  nonlinearity 
upon  a  simplified  version  of  these  equations,  which  has  the  same  stability 
properties  as  the  complete  equations.  The  results  can  then  be  extended 
qualitatively  to  the  complete  equations. 

Two  simplifications  are  made  —  the  energy  equation  is  linearized, 
and  cell  density  is  fixed  as  a  constant.  The  first  condition  has  very 
little  effect  upon  the  outcome  of  experiments,  whereas  the  second  amounts 


to  a  removal  of  the  viscosity  which  is  inherent  in  the  PIC  equations,  so 
that  the  nonlinearity  is  weaker  but  acts  in  the  same  manner.  But  these 
changes  permit  us  to  restrict  our  attention  to  a  single  nonlinear  term, 
thereby  simplifying  the  analysis  and  generalizing  the  study.  As  a  re¬ 
sult  the  conclusions  will  be  applicable  to  a  wider  class  of  nonlinear 
problems . 

Our  purpose  is  to  consider  these  equations  for  a  fluid  which  has 
been  perturbed  slightly  from  steady  state.  The  linear  part  of  the  arti¬ 
ficial  viscosity  is  set  to  zero  (so  the  system  is  unstable)  but  the  non¬ 
linear  part  is  present.  With  the  foregoing  changes  and  a  polytropic 
equation  of  state,  Eqs.  (1)  become 
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where  1^  is  the  initial  value  of  specific  internal  energy. 


Computer  Results 

A  series  of  problems  were  run  on  the  computer  in  order  to  test  the 
effect  of  the  various  parameters  on  the  final  equilibrium  attained  by 
these  equations.  An  initial  perturbation  was  supplied  as  a  small  velocity 


-19- 


in  each  cell,  and  the  resulting  fluctuations  were  observed  through  their 
effect  on  the  kinetic  energy  histories.  In  addition,  the  velocity  pro¬ 
files  were  analyzed  in  detail  through  Fourier  decomposition. 

Figure  2,  curves  (a)  to  (h),  shows  the  kinetic  energy  histories 
for  a  series  of  machine  runs  in  which  the  size  of  the  system  varies  from 
6  to  80  cells,  all  other  parameters  being  fixed  at  the  following  values: 

5t  =  0.25,  6x  =  1.0,  f  =  1.0,  7  =  1.67,  and  IQ  =  0.9.  Boundary  conditions 
specify  u  =  0  at  the  ends  of  the  region,  and  when  the  linear  form  of  vis¬ 
cosity  is  included  in  Eqs.  (7)  the  final  equilibrium  condition  is  u  =  0 

everywhere.  Thus  the  kinetic  energy  ~  is  a  good  measure  of  the 

J 

nonvanishing  fluctuation.  Since  the  frequency  of  oscillations  of  kinetic 
energy  is  too  high  to  plot  on  the  scale  of  Fig.  (2),  the  curves  trace  the 
loci  of  the  maximum  and  minimum  points.  The  average  midpoint  between 
these  lines  at  late  times  is  taken  as  the  equilibrium  value. 

In  these  curves,  notice  that  fluctuations  grow  rapidly  at  first  (as 
in  a  linear  problem)  but  are  eventually  bounded  to  equilibrium.  The  amount 
of  overshoot  which  precedes  equilibrium,  as  well  as  the  mean  equilibrium 
amplitude,  generally  increase  with  the  size  of  the  system,  although  the 
6  cell  system  is  anomalous  in  both  respects. 

The  machine  results  show  that  the  course  of  any  one  calculation 
usually  can  be  divided  into  three  phases.  In  the  first  or  linear  phase, 
the  fluctuations  grow  rapidly  in  time.  This  phase  is  terminated  by  the 
achievement,  in  the  mean,  of  velocity  fluctuations  large  enough  to  give 
appreciable  nonlinear  dissipation.  In  the  second  phase,  there  is  a  first 


TIME 

Fig.  2:  Kinetic  energy  histories  for  PIC-like  systems  which  have  been 
perturbed  from  steady  state.  The  lines  trace  the  loci  of 
maximum  and  minimum  values  of  kinetic  energy.  Input  data, 
other  than  the  number  of  cells  in  the  system,  is  the  same  for 
nil  these  problems:  &x  =  1 .0,  6t  =  0.25,  7  =  5/3,  ac^  =  0, 

f  =  1 .0,  u.  =  0.01,  I.  =  0.9. 

'  J  J 
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order  balance  between  the  instability  and  the  dissipation,  but  also  there 
is  a  higher  order  imbalance  leading  to  slow  transition  to  final  steady 
state  (the  third  phase).  The  origin  of  these  phases  will  be  discussed  in 
more  detail  in  the  following  sections.  Most  of  the  calculations  were  not 
run  long  enough  to  show  the  Phase  III  behavior. 

Other  experiments  were  performed  in  which  the  time  interval  and  the 

coefficient  of  the  nonlinear  term  were  varied.  It  was  found  that  the 

2 

final  equilibrium  kinetic  energy  was  proportional  to  (6t/f)  ,  which  be¬ 
havior  is  explained  within  the  next  few  pages. 

The  initial  perturbation  in  these  problems  was  such  that  only  the 
symmetric  (odd)  Fourier  modes  of  oscillation  were  originally  present  in 
the  system.  Furthermore,  the  rate  of  growth  of  the  even  modes  was  ob¬ 
served  to  be  very  small  compared  to  that  of  the  odd  modes,  so  that  the 
nonsymmetric  modes  never  contributed  significantly  to  the  energy  of  the 
system.  A  Fourier  analysis  of  velocity  profiles  showed  that,  among  the 
odd  modes  at  late  times,  one  mode  of  oscillation  usually  was  dominant. 
When  the  number  of  cells  in  the  system  was  small,  the  lowest  frequency 
mode  was  predominant,  exceeding  the  other  modes  in  amplitude  by  several 
orders  of  magnitude;  but  in  larger  systems  one  of  the  higher  frequency 
modes  generally  dominated.  In  these  latter  cases  the  amplitude  of  the 
dominant  mode  was  usually  two  or  three  times  that  of  the  next  largest 
mode.  Also,  although  the  wave  number  of  the  dominant  mode  differed  from 
problem  to  problem,  it  was  found  that  its  late  time  amplitude  was  always 
nearly  the  same.  These  features  are  illustrated  in  Table  1 . 


TABLE  1 


AMPLITUDE  OF  SIGNIFICANT  MODES  AT  COMPLETION 
OF  MACHINE  RUNS  FOR  SEVERAL  PROBLEMS 

The  amplitudes  shown  in  parentheses  are  averages 
of  highly  oscillating  values. 


Mode  Number  of  Cells  in  System 

No.  _ (Duration  of  Problem  in  Time  Cycles) 


12 

30 

40 

65 

80 

(9200) 

(9150) 

(13,700) 

(9700) 

(9900) 

1 

0.1  4o 

— 

0.015 

— 

— 

3 

(0.002) 

0.023 

0.014 

0.006 

— 

5 

0.000 

0.073 

0.051 

0.013 

— 

7 

0.000 

0.107 

0.121 

0.022 

— 

9 

(0.001 ) 

0.024 

0.01 1 

0.032 

0.007 

n 

(0.001) 

(0.009) 

0.001 

0.048 

0.028 

13 

(0.002) 

o.oo4 

0.108 

0.109 

15 

(o.oo4) 

0.001 

0.033 

0.032 

17 

(0.006) 

(0.002) 

0.015 

0.033 

19 

(0.010) 

(0.002) 

0.032 

0.036 

21 

(0.011) 

(0.003) 

(0.003) 

(0.009) 

23 

(0.007) 

(0.001 ) 

(0.002) 

(0.007) 

25 

(0.005) 

( 0.001  ) 

(0.002) 

(0.006) 

27 

(0.003) 

(0.002) 

(0.002) 

— 

29 

(0.002) 

0.000 

(0.003) 

— 
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The  Mechanics  of  Dissipation 


In  the  analysis  of  the  problem  of  equilibrium  balance,  it  will  turn 
out  to  be  sufficiently  accurate  to  consider  an  approximation  to  Eqs.  (7)* 
We  therefore  expand  Eqs.  (7)  in  Taylor  series  about  the  center  of  the  jth 
cell  and  about  time  n&t.  Neglecting  terms  higher  than  the  first  in  Sx  and 
5t,  we  have 


0  ox 


Now,  to  zero  order 


-(r  -  D 


so  that 


=  -(r  .  1) 


-  <r  -  Da Io^lZa. 

0  3x 


With  this  substitution  Eqs.  (8)  become 


^-1/2  _ 
“St 


=  <7  -  1}  +  h  |Lf  5x|Uj-l/2|  -  (7  -  1)2l0  1 


2  6tl  ^.1-1/2 


0  dx 


The  stability  of  Eqs.  (9)  depends  upon  the  sign  of  the  quantity 
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Tj-i/aJ "  f  5x|Vi/2 


,  ,>2.  6t 

-(r  -  i )  i0  -• 


A  stability  analysis  of  these  equations  for  a  fixed  value  of  a  will  verify 
this.  In  Eqs.  (9)  let 


Then 


n  .  ikx  cot 

"j-l/2  =  A  6  6 


n  „  ikx  cut 

1-1/2-®  C  6  • 


A[cd  +  ok  ]  +  B[ik(7  -  1)]  =  0, 
A[ik(y  -  1)L]  +  Ba>  =  0. 


For  a  nontrivial  solution  we  must  have 


2  2  2 

co(cd  +  ok  )  +  k  (7  -  1)  1^=0, 


or 


i  772 

~~  o  » 


CD  =  -  2  ok  _ 


o  k  -  tof(7  -  1) %  . 


Now  for  stability,  the  real  part  of  cd  must  he  negative.  This  will  he  true 
only  when  a  >  0. 

Since  at  steady  state 
is  initially  unstable.  But,  as  the  velocities  grow,  a  increases  and  equi¬ 
librium  is  attained  at  o  =  0.  The  corresponding  mean  velocity  magnitude 


o^u1)  definitely  negative,  the  system 


at  equilibrium  is  therefore  given  by 


u 


(7  -  1)  I0  5t 
2f  5x 


(10) 
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If  we  assume  a  probability  distribution  P(u)  of  velocities  about  this 
mean,  the  associated  mean  kinetic  energy  for  a  system  of  N  cells  of  mass 
m  each  is  given  by 


_  1  /  2 

KE  =  —  Nm  J  P(u)  u  du. 

—  CO 


(11) 


For  a  normal  distribution 

P(u)  =  0.538  o  e'°*91  U  u  , 


p  = 


f  5x 


( 7  ~  1)  IQ 


this  gives 


Nm 


KE  =  0.275  ^  =  0.275  Nm 


p 


(7-1)  IQ 

f  ox  . 


(12) 


In  contrast,  for  a  sharp  distribution  in  which  P(u)  =  &(u  -  u),  we  get 


—  1  „ 

KE  =  x  Nm 


(7  -  1)%  &t-2 


2  L  2f  6x 


=  0.125  Nm 


-(7  -  i)2i0  st^2 


f  5x 


(13) 


A  comparison  of  these  predictions  with  computer  results  for  typical 
values  of  the  parameters  is  shown  in  Fig.  3»  The  points  in  the  figure 
indicate  observed  values,  while  the  upper  line  is  a  plot  of  relation  (12), 
referring  to  normal  distribution  of  velocities  about  the  mean,  and  the 
lowest  line  shows  the  kinetic  energy.  Eg.  (13),  which  would  be  attained 
if  all  cells  had  the  mean  velocity.  It  is  seen  that  the  actual  distribu¬ 
tion  lies  between  these  extremes.  The  prediction  requires  further  analy¬ 
sis,  the  results  of  which  are  shown  by  the  middle  line  in  the  figure  and 
are  derived  below. 
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Modal  Exchange  of  Energy 

Thus  far,  the  investigation,  while  shedding  some  light  on  the  mechan¬ 
ics  of  the  damping  process  and  providing  an  estimate  of  the  equilibrium 
amplitude,  has  told  nothing  about  the  normal  modes  of  oscillation  in  the 
system  or  of  the  energy  sharing  between  these  modes.  For  this  purpose 
we  need  a  more  detailed  study  of  the  difference  equations. 

Consider  again  Eqs.  (8), 


Notice  that  we  have  dropped  the  subscripts  and  superscripts  for  ease  of 
writing. 

Following  Kryloff  and  Bogoliuboff  [3],  we  will  assume  that  the  solu¬ 
tion  of  these  first  order  equations  will  not  differ  much  from  the  solu¬ 
tion  of  the  zero  order  equations,  and  we  will  account  for  this  difference 
by  allowing  the  amplitude  and  phase  of  the  zero  order  solution  to  vary 
with  time.  An  appropriate  solution  of  the  zero  order  equations  is 

u  =  A  sin  kx  sin  (o>t  +  9), 

cos  kx  cos  (o^t  +  cp)  j 

•where 

O)  ~  ( 7  -  1  k* 

The  first  boundary  condition 
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u  =  0  at  x  =  0 

is  already  satisfied.  Applying  the  boundary  condition 


u  =  0 


at  x  =  L  = 


determines  the  unique  values 


,  _  nit  „  _  ,  0 

k  —  l  y  n  —  1  ,2, 


N  5x 

which  k  may  assume, 


Thus  the  complete  solution  of  the  zero  order  equations  can  be  written 


00 


njtx 


A  sin  — —  sin  Q  , 
.  n  L  n' 

n=1 


(H) 


1  *  ^0  An  005  TT  C0S  V 

n=1 


where 


mt  ( y  -  1  )s/l7  t 

Q  =  - - -  +  CD  =  CD  t  +  Cp  . 

hi  L  rn  n  Tn 


Now,  considering  A^  and  cpn  as  functions  of  time,  we  have 


du 

3t 


.  ,  mocv 

sin  (— ) 


L  L  n 


n«(7  -  1>/ll 

A  sin  Q  +  A  - = -  cos  Q  +  A  cp  cos  Q. 

n  n  n  n  n  n  n. 


n 

so  that,  to  impose  the  zero  order  solution,  it  is  required  that 

o  e 

A  sin  Q  +  A  9  cos  Q  =  0.  ( 1J?) 

n  n  n  n  n 

Substitution  of  solution  (14),  subject  to  condition  ( 1 5) ,  into  Eqs.  (8) 
gives 
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not (7  -  1)s/l^ 


I 


n 


nrt(7  -  1)n/57  5t 


2L 


_ _  A 

.cos  Q  n 
n 


sin  Q 


n 


sin 


notx 


La 

.  notx  .  -  1 

sm  — —  sm  Q 

E  A 

mic 

motx  .  -  \ 

cos  — —  sm  Q 

.  n  n 

L  n 

m  m 

Li 

L  my 

plus  an  identity  equation.  Thus 


mt(7  -  iy/i0  5t  T  A 


nrt(7  -  1)s/l7 

— - A  - 7 - -  sin  Q 

Lcos  Q  n  L  nj 


n 


=  f  Bx  /  sin 
J0 


B=  A  (  E  A  sin  is  sin  0 
L  ox  \  i  |  L  J 


mat  matx 

.  —  cos  — — 
m  I  L 


\ 


V  a  mat  matx  .  „ 

LA  —  cos  — —  sm  Q  J 


Integration  of  the  right-hand  side  hy  parts  and  simplification  give  an  ex¬ 
pression  for  the  rate  of  growth  of  the  amplitudes, 

1 


/not(7  -  1  y/T\  4f 

A  =  A  - 3= - 2  J  sin  2Q - - 

n  nV  2L  y  n  L(7-1)^5t 


cos  Q 


/ 

nJQ 


cos 


notx 


.  .  iotx  .  _ 

A  sm  — —  sm  Q 

*  Xj  / 


\ 


v  a  mjt  mjtx  .  .  \  - 

LA  —  cos  — —  sm  Q  )  dx. 
^  m  L  L  W 


(16) 


Consider,  first,  the  case  f  =  0.  Combining  Eqs.  (15)  and  (l 6)  we 
obtain 

cp  =  -  o)  sin  (co  t  +  cp  ) , 
rn  n  n  n 

which  on  integration  gives 

cp  =  -  cn  t  4*  tan”1  (a)t  +  C  ),  ( 1 7) 

vn  n  n  n  7 

where  is  an  arbitrary  constant.  With  this,  then 


dx. 
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on  t  +  C 
n  n 


A 


n 


cd  A 
n  n 


L1  + 


(cd  t 
'  n 


C  )2j 
n 


Thus 


A 

n 


+  (co  t  +  c  )  , 

n  n 


(18) 


where  K  is  a  second  arbitrary  constant, 
n 

This  shows  that  with  f  =  0  ve  can  expect  a  growth  in  amplitude  for 
each  component,  becoming  linear  in  time  for  large  times.  This  solution 
is  also  appropriate  for  f  /  0  when  the  amplitudes  are  small  enough  to 
neglect  the  nonlinear  term  in  Eq.  (16).  But  as  the  amplitudes  grow,  this 
nonlinear  term  will  eventually  check  the  instability,  bringing  the  system 
to  equilibrium. 

Consider  now  the  case  f  /  0.  By  neglecting  cross  product  terms  (whose 
contribution  is  small  for  the  significant  lower  frequency  modes),  we  can 
write  Eq.  (16)  in  a  simpler  form  which  illustrates  the  manner  in  which 
fluctuations  are  damped. 


A  =  A  cd  sin  Q  cos  Q  1 
n  n  n  n  n  L 


4f  5x _ 

L(r  -  i)Vi^  st 


2 

cos 


njtx 

L 


(16*) 


The  amplitudes  increase  in  magnitude  until  the  velocity  becomes  large 
enough  to  make  the  bracket  term  small.  Thus,  the  mean  value  of  A^  ap¬ 
proaches  zero  to  first  order  and  An  achieves  its  maximum  value.  The  order 
in  which  the  modes  are  maximized  depends  upon  the  configuration  of  the  ve¬ 
locity  profile  and  of  this  we  can  say  very  little  at  first.  But  the  ve¬ 
locity  magnitude  must  continue  its  growth  as  long  as  there  is  a  single 
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O 

inode  for  which  A  ^  0}  during  "this  period  "the  bracket  "berms  become  nega- 
n 

tive  for  many  modes  and  these  oscillations  decay. 

Herein,  perhaps,  lies  the  explanation  of  the  origin  of  the  dominant 
mode  oscillation,  which  was  discussed  earlier.  For  consider  the  system 
at  the  time  when  there  is  but  a  single  mode  which  remains  to  be  maximized. 
The  velocity  is  increasing  in  magnitude  but  all  other  frequency  oscilla¬ 
tions  are  declining,  and  hence  the  velocity  profile  is  approaching  closer 
and  closer  to  the  configuration  of  the  growing  oscillation.  It  can  be 
shown  that  the  integral  in  Eq.  (l 6*)  is  least  when  the  velocity  is  com- 

Q  ___ 

posed  entirely  of  — —  frequency  oscillations;  therefore  the  growing  re¬ 
semblance  of  the  velocity  to  this  frequency  only  serves  to  prolong  the 
growth  of  this  final  mode  and  the  decay  of  all  other  modes.  Thus  the  end 
of  the  first  phase  of  the  process,  which  corresponds  to  the  time  of  maxi¬ 
mum  velocity  magnitude,  finds  a  major  portion  of  the  energy  of  the  system 
concentrated  into  a  single  mode;  this  concentration  increases  throughout 
Phase  II. 

To  include  the  cross  product  terms  and  perform  the  analysis  in  general 
would  be  difficult;  however,  when  a  particular  mode,  say  number  a,  domi¬ 
nates  the  system  to  the  extent  that 

X  ki sin  nr sin 

l 

then  considerable  simplification  is  possible*  Making  use  of  this  "domi¬ 
nant-mode n  assumption,  which  will  he  done  for  the  rest  of  this  paper, 

Eq.  ( 1 6)  can  he  written 


A 

a 


sm 


a  jtx 


sm 


(19) 
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In  the  machine  calculations,  it  was  observed  that  the  even  modes 
never  contributed  significantly  to  the  fluctuation  energy.  This  is  rea¬ 
sonable  to  expect,  since  the  equations  with  f  =  0  conserve  symmetry  and 
initially  only  odd  modes  were  present;  by  the  time  even  modes  could 
couple  in  significantly,  equilibrium  of  Phase  II  type  had  been  achieved 
and  further  even-mode  growth  was  slowed  to  second  order.  Thus,  with  only 
odd  modes  considered,  the  sum  in  Eq.  (20)  can  be  written 
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> 

n  Z_j  m 


sin 


a  irx 


m 


*1 


mjtx  nj rx  , 

cos  — —  cos  — —  ax 


2a  -  n|  v2 a  +  n  ^|4a  -  n|  ¥4a  +  n  _ 


'V! 


L  , 

”  n  Kfn  3  3  15  "  15 


,  ligag  -  nl  ,  llgg£L±-5i  ...\ 

+  p  +  2  /• 

1  -  4q  1  -  4q 


When  n  =  a,  the  dominant  mode,  only  the  first  two  terms  in  the  sum  contri¬ 
bute  significantly  to  S  .  Thus 


S  -»•§  —  aA  sin  Q  . 

a  3  it  a  3  a  a 

When  n  4  a.  we  will  assume  that  S  can  be  limited  to  its  first  term.  This 
assumption  will  be  justified  presently,  at  least  for  the  significant  lower- 
number  modes,  i.e.,  for  n  <  2a.  With  this  assumption  the  system  of 
Eqs.  (20)  can  be  written 


A  =  A  a) 

a  a  a 


A  03 

n  n 


sin  Q  cos  Q 

a  a 

sin  Q  cos  Q 
n  n 


(’ 

0 


2 

3 


4f  5x _ _ 

*(7  -  1 )\  &t 

4f  5x _ 

* ( 7  “  1 ) 2I0  St 


a 


a 


|3in 

lBln  tl)' 


(21) 

n  j-  a. 


Since  these  equations  have  been  derived  on  the  basis  of  the  dominant 
mode  assumption,  which  became  effective  at  the  close  of  Phase  I  of  the 
process,  they  are  appropriate  for  a  description  of  quasi-equilibrium  and 
equilibrium  conditions,  that  is.  Phases  II  and  III.  Without  solving  them 
we  may  notice  their  late-time  properties.  The  rate  of  change  of  all  the 
modes  depends  upon  the  magnitude  of  A  which  will  continue  to  increase. 
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on  the  average,  as  long  as  the  bracket  term  in  the  first  equation  is  posi¬ 
tive.  Before  this  growth  stops,  however,  the  bracket  term  in  the  equation 


for  A  .  n  4  a.  will  have  become  negative.  Thus  the  final, 
n  ' 

librium  should  correspond  to 

3  3*2(7  -  1)2IQ  6t 

Aa  2  4f  5x  <  |  sin  |  >  "  16f  6x  * 

An  0#  n  4  a. 


Phase  III,  equi 


(22) 


where  <  >  signifies  an  average  over  time.  However,  there  is  some  evidence 
to  indicate  that  the  equilibrium  amplitude  for  may  be  changed  somewhat 
by  the  effect  of  the  sinusoidal  terms  in  the  expression  for  A^.  This  mat¬ 
ter  will  be  returned  to  again  below. 

Phase  II  is  that  period,  following  the  initial  amplitude  growth,  when 
the  system  is  approaching  this  single  frequency  oscillation.  Its  duration 
varies  considerably  from  problem  to  problem  and  is,  as  shall  be  demon¬ 
strated,  largely  dependent  upon  the  number  of  the  dominant  mode.  Further¬ 
more,  the  number  of  the  dominant  mode  increases  with  system  size,  so  that 
ultimately  the  rate  of  energy  concentration  depends  upon  the  system  size. 
This  fact  is  very  much  in  evidence  in  the  kinetic  energy  profiles  of 
Fig.  2.  The  6,  12,  and  18  cell  systems,  for  which  Fourier  decomposi¬ 
tions  indicate  the  lowest  frequency  mode  is  dominant,  pass  rapidly  from 
initial  damping  to  the  Phase  III  equilibrium  condition  of  a  uniform  ampli¬ 
tude,  single  frequency  oscillation.  But  the  larger  systems,  which  are 
dominated  by  higher  frequency  modes,  exhibit  profiles  indicative  of  com¬ 
posite  frequency  oscillations.  Furthermore,  the  amplitude  of  these 


kinetic  energy  curves  increases  with  time,  indicating  a  slow  growth  in 
the  magnitude  of  the  dominant  mode  at  the  expense  of  the  secondary  modes. 

The  variation  in  the  calculation  time  required  to  attain  a  true 
dominant  mode  distribution  of  energy  and  the  dependence  of  this  calcula¬ 
tion  time  upon  system  size  are  apparent  from  Table  1 ,  which  shows  the 
amplitude  of  various  modes  at  the  completion  of  the  machine  runs.  This 
completion  time  is  quite  arbitrary  in  that,  although  all  the  problems 
have  passed  the  stage  of  initial  amplitude  growth,  there  is  a  marked  con¬ 
trast  in  the  amount  of  energy  concentration  which  has  taken  place.  The 
12  cell  system  is  the  only  one  in  which  essentially  all  the  energy  of  the 
system  has  been  concentrated  in  one  mode.  The  other  problems,  one  of 
which  was  run  considerably  longer  than  the  12  cell  problem,  all  contain 
secondary  modes  of  significant  amplitude. 

We  have  not  been  able  to  predict  a  priori  which  mode  will  dominate 
nor  to  explain  why  the  number  of  the  dominant  mode  increases  with  the  num¬ 
ber  of  cells  in  the  system;  but  we  can  show  qualitatively  why  the  secon¬ 
dary  mode  amplitudes  recede  more  slowly  when  the  dominant  mode  wave  number 
increases.  For  this  purpose  it  seems  appropriate  to  make  use  of  Eqs.  (21), 
even  though  the  system  has  not  yet  reached  a  true  dominant  mode  condition. 
They  are  considered  applicable  because  the  amplitude  of  the  dominant  mode 
is  so  much  larger  than  the  secondary  mode  amplitudes  that  condition  (19) 
should  be  valid.  Simplifying  Eqs.  (21)  we  write 

A  =  A  o>  (l-fgAj, 

a  a  a  \  3  as 


A  =  A  a)  l  1  -  gA  J, 
n  n  n  v  'v// 


where 


g  = 


4f  6x  <  | sin  Qa  | > 

*(7-1)%  st 


Solving  the  first  equation,  we  obtain 
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where  R  is  an  arbitrary  constant.  Substituting  this  expression  for  A 


a 


into  the  second  equations  yields  a  solution  for  Aq, 


A  =  K 
n  n 


-.n/a 


CD 


a 


5/2 


(,  ♦  |  k  Be  a 


(23) 


where  K  is  another  arbitrary  constant, 
n 


Notice  that  as  t  -*  °°, 


Aa  2g' 


n 


0, 


so  that  the  final  equilibrium  solution  is  the  same  as  before.  But  Eq.  (23) 
indicates  why  this  final  solution  is  delayed  as  the  value  of  a  increases. 
When  a  =  1,  as  in  the  12  cell  problem,  the  smallest  possible  value  of  the 
exponent  is  three,  so  that  all  the  secondary  modes  lose  amplitude  rapidly; 
but  for  larger  values  of  a  there  exist  modes  for  which  this  exponent  is 
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approximately  one,  so  that  the  decay  of  these  amplitude  curves  is  much 
more  gradual.. 

Likewise  this  equation  demonstrates  why,  for  a  fixed  value  of  a, 
the  high  frequency  vibrations  are  damped  much  more  rapidly  than  the  lower 
frequency  ones,  substantiating  our  assumption  that  the  high  frequency  con¬ 
tributions  to  Sn  could  be  neglected  when  n  <  2a .  This  rapid  decay  of  the 
high  frequency  modes  is  apparent  in  Table  2,  which,  for  a  system  of  40 
cells,  compares  the  amplitude  of  significant  modes  at  an  earlier  Phase  II 
time  with  those  at  problem  completion  time  as  listed  in  Table  1 .  Notice 
also,  in  Table  2,  the  growth  of  the  dominant  mode  amplitude. 


TABLE  2 

AMPLITUDE  OF  SIGNIFICANT  MODES  AT  EARLY  AND  LATER 
EQUILIBRIUM  TIMES  FOR  A  SYSTEM  OF  40  CELLS 


Time  Mode  Number 


1 

3 

5 

7 

9 

11 

13 

7,300 

0.015 

0.016 

O.053 

0.093 

0.061 

0.022 

0.017 

13,700 

0.015 

0.014 

0.051 

0.121 

0.011 

0.001 

o.oo4 

Now  as  the  amplitude  of  these  high  frequency  vibrations  begins  to 
change  rapidly,  their  phase  angles  once  again  assume  a  strong  time  depend¬ 
ency  as  a  result  of  the  condition  expressed  in  Eq.  ( 1 5) *  These  complica¬ 
tions  are  made  apparent  in  the  computer  results  by  rather  erratic  varia¬ 
tions  in  the  amplitude  and  period  of  these  high  frequency  oscillations  at 
late  times.  The  irregular  character  of  these  oscillations  is  indicated  in 
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Table  1  by  parentheses  around  their  average  amplitude  value. 

On  the  other  hand,  in  the  larger  systems,  the  lower  frequency  oscil¬ 
lations  are  extremely  uniform  in  both  amplitude  and  period  at  the  time 
that  these  computer  runs  were  completed.  This  would  indicate  that  the 
dominant  mode  amplitude  has  not  as  yet  changed  much  from  the  value  which 
would  make  =  0  in  Eq.  (21),  i.e., 


Ay  -  i)2i0  6t 

Aa  8f  5x 


(24) 


And,  indeed,  the  dominant  amplitudes  in  Table  1  do  not  vary  much  from 
this  value;  the  greatest  variation  is  13$* 

Since  this  rather  slow  rate  of  growth  of  the  dominant  mode  is  some¬ 
what  at  odds  with  what  one  would  predict  from  Eqs.  (21),  it  is  perhaps 
worthy  of  some  additional  comment.  Notice  in  these  equations  that  the 

o 

expression  for  A  is  a  product  of  two  factors. 


ka  “a  sin  %  oos 


and 


0 


2  4f  5x 
5  n(y  -  1)2Iq  Bt 


Aa  |sin  Qa|J  • 


The  term  in  parentheses  is  positive  for  the  Phase  II  dominant  mode  ampli¬ 
tudes  given  by  Eq.  (24),  and  hence  any  slowing  down  of  the  rate  of  growth 
must  result  from  the  first  term.  This  could  happen  if  the  first  term  was 
approaching  an  oscillating  function  of  time,  in  which  case  the  product  of 
the  two  terms  would  oscillate  and  would  experience  alternating  periods 
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of  growth  and  decay.  As  this  first  term  becomes  a  sinusoidal  function, 

o 

A  vanishes  on  the  average  and  the  system  reaches  final  equilibrium o 
a 

There  is  some  reason  to  believe  that  the  final  equilibrium  state, 
which  these  machine  problems  are  approaching,  is  characterized  more  by 
the  sinusoidal  nature  of  this  first  term  than  by  the  vanishing  of  the 
term  in  parentheses.  Evidence  in  support  of  this  is  given  by  the  ampli¬ 
tude  of  the  dominant  mode  of  the  12  cell  system  at  problem  completion 
time.  Both  the  extreme  concentration  of  energy  into  a  single  mode,  evi¬ 
dent  in  Table  1 ,  and  the  uniform  kinetic  energy  amplitude  of  Fig.  2b  in¬ 
dicate  that  this  12  cell  problem  is  at  least  very  close  to  a  final  equi¬ 
librium  state.  And  yet  the  dominant  mode  amplitude  for  this  problem  is 
0.140,  which  is  only  7 6$  of  the  value  that  would  make  the  term  in  paren¬ 
theses  in  Eq.  (21)  vanish.  Hence  the  constancy  of  the  dominant  mode 
amplitude  must  result  from  the  vanishing  on  the  average  of  the  first  term. 

This  matter  will  be  referred  to  again  in  the  discussion  of  final 
kinetic  energy  predictions. 


Predictions 

Consider,  now,  an  estimate  of  the  total  kinetic  energy  of  the  system 
on  the  basis  of  this  Phase  II  dominant  mode  amplitude  [Eq.  (24)].  The 
kinetic  energy  in  the  jth  cell  is  given  by 
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so  that,  with  m  =  1 , 
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This  prediction  is  shown  as  the  middle  line  in  Fig.  3*  The  agreement  is 
considerably  better  (at  this  stage)  than  that  obtained  in  the  first  analy¬ 
sis;  most  of  the  discrepancy  can  be  attributed  to  the  actual  strength  of 
the  neglected  secondary  modes. 

But  from  this  intermediate  stage  of  equilibrium  we  expect  the  kinetic 
energy  profile  to  rise  rather  slowly  as  the  dominant  mode  amplitude  ap¬ 
proaches  its  asymptotic  value  given  by  Eq.  (22).  Since  the  asymptotic  dom¬ 
inant  mode  amplitude  is  3/2  the  value  used  in  the  kinetic  energy  determina¬ 
tion  above,  we  might  expect  that  the  kinetic  energy  of  the  system  would 
eventually  attain  a  level  9/4  higher  than  this  intermediate  plane.  How¬ 
ever,  it  is  probable  that  this  growth  will  be  tempered  by  the  fact  that, 

•  Q 

as  A  becomes  small,  the  first  factor  in  the  expression  for  A  will  become 

ot  ^ 

sinusoidal,  so  that  will  become  an  oscillating  function  of  time  in 
Eq.  (21).  The  final  kinetic  energy  level  should  then  lie  at  some  inter¬ 
mediate  plateau,  probably  not  far  removed  from  the  prediction  made  by  the 
first  analysis.  Unfortunately,  a  prohibitive  amount  of  machine  time  would 
be  required  to  enable  one  of  these  large  systems  to  reach  this  ultimate 
goal. 
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It  is  also  possible  to  obtain  a  fairly  accurate  estimate  of  the  fre¬ 
quency  of  the  dominant  mode  at  equilibrium.  When  the  system  has  attained 
a  true  dominant  mode  energy  distribution,  indicative  of  Phase  III  equili¬ 
brium,  and  the  dominant  mode  is  approaching  its  limiting  value,  then  its 

frequency  should  be  approaching  the  natural  frequency,  since  cpa  -» 0  with 

• 

Aq,  according  to  Eq.  (15).  In  the  12  cell  system,  which  is  a  case  of 
this  type,  the  frequency  of  the  first  mode  differs  by  less  than  2$  from 
its  natural  frequency. 

But  for  a  system  which  has  only  attained  the  intermediate  stage  of 
equilibrium,  a  prediction  is  somewhat  more  difficult  to  obtain.  The  re¬ 
lationship  between  the  amplitude  and  phase  of  the  dominant  mode  in  such 
a  system  is  described  by  Eqs.  (21)  and  (15)# 
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sin  Q 
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8a  f  5x _ 
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A  sin  Q  +  A  (Q  -  cd  )  cos  Q  =  0. 
a  a  a  a  a  a 

Since  the  experimental  evidence  indicates  that  the  amplitude  and  period 
of  the  dominant  mode  vary  quite  slowly  at  this  stage,  let  us  set 


a 


K  +  Z, 


=  fit  +  e 


in  these  equations,  where  £  and  e  are  higher  order  correction  terms. 
Neglecting  higher  order  terms,  we  get 

£  =  ^  K  sin  2fit  (o)Q  -  ZK  |  sin  fit  | ) , 


-k2- 


I  sin  Qt  +  K(fl  -  flcu,  +  e)  cos  (at  +  e)  =  0, 


a 


or 


sin  2ftt  (cua  -  ZKjsin  flt|)  +  a  -  (oq  +  €  =  0. 


Now  average  over  time  and  assume  that  <  e  >  -  0  to  get 
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For  K  use  the  Phase  II  dominant  mode  amplitude  given  by  Eq«  (24).  Then 
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Table  3  shows  a  comparison  of  this  predicted  dominant  mode  frequency  with 
the  observed  results  for  problems  at  this  intermediate  stage  of 
equilibrium . 

TABLE  3 


COMPARISON  BETWEEN  OBSERVED  AND  PREDICTED  DOMINANT  MODE 
FREQUENCIES  AT  THE  INTERMEDIATE  STAGE  OF  EQUILIBRIUM 


Frequency 

Number  of 

Cells 

30 

40 

65 

80 

Observed 

0.427 

0.331 

0.370 

— 

to 

II 

RS 

0.439 

0.329 

0.375 

0.323 
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Application  to  the  PIC  Equations 

The  results  which  have  been  obtained  up  to  now  have  been  applicable 
to  a  simplified  version  of  the  PIC  equations.  We  should  like  to  investi¬ 
gate  the  extent  to  which  these  results  apply  to  the  complete  PIC  equations. 

The  primary  simplification  made  was  the  abandoning  of  the  motion  of 
particles  entirely  and  the  assumption  that  every  cell  had  constant  density; 
the  resulting  equations  are  expressed  as  Eqs.  (l).  Then,  for  ease  of 
analysis,  these  equations  were  further  simplified  by  linearizing  the  energy 
equation. 

Let  us  consider,  first,  the  effect  of  this  latter  modification  on 
equilibrium  attainment.  Figure  4  shows  the  kinetic  energy  histories  for 
two  problems,  both  of  which  consisted  of  40  cell  systems  with  the  same 
input  data  as  the  problems  in  Fig.  2.  The  dashed  lines  in  the  figure  dup¬ 
licate  Fig.  2e,  while  the  solid  lines  represent  the  maximum  and  minimum 
kinetic  energy  values  for  a  problem  in  which  the  nonlinear  terms  in  the 
energy  equation  were  retained.  The  average  late  time  equilibrium  kinetic 
energy  is  O.O83  for  the  linearized  version  as  opposed  to  0.081  for  the 
nonlinear  form.  The  reason  why  the  discrepancy  is  so  slight  is  that  each 
term  in  the  nonlinear  energy  equation  is  of  higher  order  in  the  perturba¬ 
tion  than  the  corresponding  term  in  the  momentum  equation.  Because  the 
perturbation  is  bounded  to  a  very  low  level,  the  effect  of  these  higher 
order  terms  is  negligible. 

But  the  effect  of  the  other  modification,  the  removal  of  particles, 
was  much  more  pronounced.  For,  not  only  was  the  inherent  viscosity  of 
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the  method,  which  results  from  particles  crossing  cell  boundaries,  lost 
to  the  system  but,  in  addition,  the  system  was  deprived  of  the  damping 
effects  of  momentum  and  energy  conservation.  The  result  is  that  the 
energy  profiles  of  a  PIC  problem  are  bounded  to  a  much  lower  level  than 
the  foregoing  analysis  would  indicate  and,  furthermore,  the  damping 
effect  is  apparent  much  sooner  than  was  the  case  for  the  simplified 
method.  These  improvements  are  very  much  in  evidence  in  Fig.  5,  which, 
shows  the  kinetic  energy  profiles  for  a  PIC  problem  of  40  cells  with 
the  same  input  as  the  problems  of  Fig.  4  and  with  4  particles  per  cell. 

But  the  presence  of  particles  is  not  enough  to  insure  effective 
damping  of  fluctuations;  the  number  of  particles  per  cell  is  also  ex¬ 
tremely  important.  There  seems  to  be  an  optimum  number  of  particles 
per  cell,  below  which  effective  damping  is  not  attained  and  above  which 
only  minor  improvement  is  observed.  In  illustration.  Table  4  lists  the 
equilibrium  kinetic  energy  at  late  times  of  2,  4,  and  8  particles  per 
cell  versions  of  the  40  cell  problem  illustrated  in  Fig.  5. 

TABLE  4 

EQUILIBRIUM  KINETIC  ENERGY  AMPLITUDES  OBSERVED  FOR  40  CELL  SYSTEMS 
CONSISTING  OF  2,  4,  and  8  PARTICLES  PER  CELL 

Particles  per  Cell  Equilibrium  K«  E» 

2  ~  0.3 

4  0.020 

8  0.014 
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Fig.  5:  Kinetic  energy  history  for  a  40  cell  PIC  system \hich  has  been  perturbed  from  steady 
state.  Each  cell  contains  k  particles  of  mass  0.25;  otherwise  input  is  the  same  as 
that  of  Fig.  2e. 


However,  it  should  not  be  thought  that  the  requirement  that  a  sys¬ 
tem  contain  a  minimum  number  of  particles  per  cell,  on  the  average,  con¬ 
stitutes  a  new  restriction  on  the  PIC  method.  It  merely  reflects  the 
need  for  proper  resolution  of  the  fluid  within  the  cells. 

Once  this  basic  requirement  is  satisfied,  it  appears  that  the  rela¬ 
tionship  between  equilibrium  kinetic  energy  and  the  parameters  f,  5t,  and 
5x,  which  was  obtained  from  the  analysis  of  the  simplified  equations 
[Eq.  (25)],  holds  at  least  qualitatively  for  the  PIC  equations  as  well. 
But,  whereas  that  analysis  predicted  that  the  equilibrium  kinetic  energy 
would  vary  quadratically  with  the  ratio  5t/(f  6x),  the  PIC  experiments  in¬ 
dicate  that  the  variation  is  slower  than  that  and  depends  upon  the  param¬ 
eter  which  is  changed  in  the  ratio.  Furthermore,  in  a  PIC  problem  the 
permissible  fluctuation  of  kinetic  energy  (and  therefore  the  size  of  the 
ratio)  is  limited  by  the  fact  that  velocities  must  remain  well  below  the 
magnitude  by  which  particles  can  be  moved  one  cell  length  in  a  time  cycle, 
i.e.,  |u|  <  5x/ 6t .  Once  this  condition  is  violated,  the  fluctuations  will 
no  longer  be  bounded  in  amplitude. 

Figure  6  shows  the  kinetic  energy  histories  for  a  series  of  problems 
in  which  the  ratio  above  is  halved  by  varying  each  of  the  parameters  in 
turn.  The  top  profile  is  a  duplicate  of  that  in  Fig.  5  except  that  the 
vertical  scale  is  doubled  in  Fig.  6.  The  next  three  profiles  in  this 
figure  demonstrate  the  effect  on  kinetic  energy  of  doubling  the  viscosity 
coefficient,  f,  and  the  cell  length,  5x,  and  of  halving  the  time  interval, 
6t,  respectively.  Note  that  the  change  in  the  time  increment  had  a  much 
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f  =2.0 


8x  =  2.0 


8t  =  O.I25 


This  series  demonstrates  the  effect  of  various  parameters  on 
the  kinetic  energy  history  of  a  40  cell  PIC  system  which  has 
been  perturbed  from  steady  state.  Except  for  the  changes 
noted,  the  input  data  is  the  same  as  that  for  Fig.  2e,  in  par 
ticular,  f  =  1.0,  6x  =  1.0,  &t  =  0.25. 


greater  effect  on.  the  kinetic  energy  profile  than  either  of  the  other 
parameter  changes. 

The  importance  of  the  time  increment  on  instability  bounding,  as 
compared  to  the  effect  of  the  other  parameters,  persists  even  to  zero 
values  of  the  artificial  viscosity.  This  is  evident  in  Fig.  7  which 
shows  the  effect  of  doubling  the  cell  length  and  halving  the  time  inter¬ 
val  in  a  system  in  which  the  only  viscosity  present  is  the  effective 
viscosity  of  the  PIC  method  (Appendix  II ) . 

A  comparison  of  Figs.  6  and  7  is  demonstrative  of  the  fact  that  the 
introduction  of  nonlinear  artificial  viscosity  into  a  PIC  calculation 
will  not  accomplish  a  great  deal  in  the  way  of  instability  bounding. 

This  means  of  controlling  fluctuations  in  what  should  be  a  stagnant  re¬ 
gion  of  the  fluid  should  especially  be  avoided  in  problems  where  one  ex¬ 
pects  large  velocity  gradients.  A  far  better  way  of  controlling  these 
fluctuations  is  to  place  an  upper  limit  on  the  time  increment.  Once  the 
other  parameters  of  the  system  have  been  chosen,  this  limit  may  be  deter¬ 
mined  by  introducing  a  perturbation  into  an  otherwise  stagnant  system  and 
varying  the  size  of  the  time  interval  until  a  tolerable  kinetic  energy 
level  is  obtained.  In  many  cases  the  limit  on  the  time  increment  thus  de> 
termined  is  less  restrictive  than  that  required  to  preserve  stability  in 
the  regions  of  high  velocity  flow. 


-50- 


APPENDIX  I.  THE  FORM  OF  THE  PIC  ENERGY  EQUATION 


In  the  original  PIC  method  the  energy  difference  equation  was  derived 
from  the  differential  equation 

SE  SE  _  S(pu)  (a  d 

pdt  +  pu^ - tT’  U  } 

where  E  is  the  total  energy  (the  other  quantities  were  defined  in  Part  I). 
Using  this  difference  scheme  the  method  was  nonconservative  of  energy  in 
Phase  I. 

Energy  conservation  can  be  obtained  in  the  PIC  method  by  using  as  a 
basis  for  the  energy  difference  equation  the  differential  equation 


Si  Si 

p5t +  pu  s 


which  is  equivalent  to  Eq.  (A-l).  Dropping  the  transport  term  in  this 
equation  and  transforming  to  difference  notation  we  have 


si,h/2  .  pj-yg  L 

5x  St  25x  Vj-V2 


where  N^  represents  the  number  of  particles  in  cell  j  and  m  is  the  mass 
of  a  particle.  The  velocity,  u,  in  this  equation  is  then  replaced  by  its 
time  average  over  Phase  I,  u.  With  this  change  the  Phase  I  momentum  and 
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energy  equations,  with  artificial  viscosity  terms  neglected,  appear  as 
follows: 


(A-2) 


~  _n  &t  /  n  n) 

Uj-l/2  “  Uj-l/2  2Knm  VPj“3/2  "  Pj+l/2> 
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~  ^n  .  6t  Pj-1/2  A- 

Ij~ V 2  Io-l/2  +  ^  VUj-3/2  "  Uj+l/2/# 

j 

where  the  tilde  symbolizes  the  fact  that  these  terms  do  not  as  yet  in¬ 
clude  the  transport  effect. 

To  see  that  this  differencing  technique,  together  with  the  proper 
boundary  conditions,  does  conserve  total  energy,  observe  that  at  time  t 
the  total  energy  of  the  system  is 

En  =  I,  ">  +  i  (uj-i/a)  ]* 


whereas  by  the  end  of  Phase  I  it  has  changed  to 
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the  total  change  being 
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after  rearrangement  of  terms.  This  equation  is  now  in  conservative  form 
and,  in  terms  of  the  fictitious  cells  beyond  the  system  (j  =  0  and 
j  =  J  +  1 ) ,  can  be  written 

X 

6t  /  n  —  n—  n—  n  —  \ 

&E  =  ~2  VPl/2U-l/2  +  P-l/2Ul/2  "  P J+1  / 2UJ- 1/2  "  P J- 1  / 2  J+l/2/* 


Since  the  system  is  considered  to  have  rigid  boundaries  the  velocity  must 
vanish  at  both  ends;  thus  u_(^2  =  -U1 /2  anb  UJ-1/2  =  "*'UJ+1 /2  a^ways* 
Likewise  du/dt  =  0  at  each  boundary  so  that  p_  ^2  =  p^  2  and  pJ_1y2 


=  pJ+l/2  by  E<1°  ^A"2^* 
ary  conditions. 


Energy  conservation  then  follows  from  these  bound- 
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APPENDIX  II.  ARTIFICIAL  VISCOSITY 


Inherent  in  the  PIC  method  is  a  diffusion  term  which  is  proportional 
to  velocity  magnitude.  This  was  demonstrated  in  Reference  1,  p.  15  ff, 
by  forming  the  Taylor  expansions  of  the  difference  equations  which  repre¬ 
sent  «l 1  three  phases  of  the  PIC  method.  The  resulting  equations,  with 
higher  order  terms  neglected  and  subscripts  dropped,  were 


where 

X'  =  ~  &x  pu 

is  called  the  effective  viscosity  coefficient.  Comparing  these  equations 
with  the  momentum  and  energy  equations  expressing  one  dimensional  compres¬ 
sible  fluid  flow, 


chi 


du 


p  +  pu  5E  = 

p  +  pu  ^ 


dp 

du 

“  Pd^ 


demonstrates  that  the  X'  terms  have  their  origin  in  the  difference  method. 
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The  X'  terms  have  their  principal  effectiveness  in  regions  where  the 
velocity  gradient  is  large  or  where  fluctuations  occur  in  high  speed 
material.  In  fact,  it  is  due  to  the  presence  of  these  terms  that  the  PIC 
method  achieves  success  in  approximating  high  velocity  fluid  flow.  This 
measure  of  success  has  not  been  obtained  in  problems  in  which  the  fluid 
speed  is  considerably  subsonic. 

Specifically,  an  important  limitation  of  the  PIC  method  exists  when 
perturbations  are  introduced  into  an  otherwise  stagnant  fluid.  Large 
fluctuations  about  the  proper  values  develop  in  the  region  of  the  per¬ 
turbation.  Damping  of  these  fluctuations  will  only  develop  as  the  ve¬ 
locity  increases  to  the  point  where  the  effective  viscosity  becomes  im¬ 
portant.  In  order  to  obtain  rapid  damping  in  a  nearly  stagnant  fluid,  a 
dissipative  force  not  proportional  to  velocity  is  required. 

Provision  has  been  made  in  the  PIC  method  to  introduce,  artifically, 
variable  amounts  of  damping  forces  in  the  form  of  pressure  modifications. 
This  dissipative  term  is  of  the  form 


n 

Pj  aC 


f 

0  +  2 


Vl/2  +  Vl/2 


n 


and  is  always  calculated  at  cell  boundaries. 

The  acQ  part  of  q  is  of  the  type  discussed  above  since  it  is  not  pro¬ 
portional  to  velocity.  It  is  similar  to  a  Landshoff  [4]  artificial  viscosity 
in  that  cQ  is  taken  as  a  representative  sound  speed  for  the  system.  This 
portion  of  the  dissipative  force  is  applied  only  when  the  system  is  exper¬ 
iencing  compression. 
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In  Part  I  it  was  demonstrated.  (p.  1 5)  that  the  PIC  equations  were  un¬ 


conditionally  unstable  for 


ac 


0 


< 


+  cos 


X  (7  -  I)2  6t  I0 
N /  25x 


This  instability  is  bounded  in  amplitude,  however,  by  the  effective  viscos¬ 
ity  of  the  PIC  method.  A  comparison  of  the  recovery  from  a  density  per¬ 
turbation  away  from  steady  state  is  shown  in  Fig.  A-1  for  systems  which 
contain  acQ  type  viscosity,  f  type  viscosity,  and  effective  viscosity  only. 
Notice  that  recovery  is  complete  in  the  presence  of  acQ  type  viscosity  but 
that  a  finite  equilibrium  amplitude  is  approached  in  the  other  two  cases 
(see  Part  II). 

The  f  part  of  the  artificial  viscosity  is  proportional  to  the  effec¬ 
tive  viscosity.  It  may  be  applied  either  in  compressions  and  rarefactions 
or  in  rarefactions  alone.  Thus  one  is  permitted  a  rather  wide  choice  in 
either  bolstering,  reducing,  or  removing  the  effective  viscosity  of  the 
system.  Such  a  capability  is  important  in  those  problems  in  which  effec¬ 
tive  viscosity  has  a  detrimental  effect.  As  an  example,  when  gas  moves 
away  from  a  wall  the  effective  viscosity  can  cause  cavitation  by  smoothing 
the  velocity  profile.  In  this  case  removing  effective  viscosity  only  in 
rarefactions  prevents  cavitation  while  retaining  the  smoothing  effect  in 
compressions. 

However,  in  most  problems  in  which  a  rather  wide  range  of  velocity 
values  is  encountered,  better  smoothing  seems  to  be  obtainable  from  the 
acQ  type  of  viscosity  than  from  any  other  type.  As  an  example,  when  a 
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Fig.  A-1 :  Profiles  of  maximum  kinetic  energy  from  three  20  cell  PIC 
systems  which  have  been  perturbed  from  steady  state  as  the 
result  of  an  instantaneous  density  imbalance.  The  dotted 
line  profile  resulted  from  calculations  employing  linear 
artificial  viscosity  (ac^  =  1.0),  the  dashed  line  profile 
corresponds  to  nonlinear  artificial  viscosity  (f  =  1.0), 
and  the  solid  line  to  a  system  which  did  not  employ  arti¬ 
ficial  viscosity. 
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plane  steady  shock  strikes  and  reflects  from  a  rigid  vail,  one  would  ex¬ 
pect  the  material  velocity  in  the  vicinity  of  the  wall  to  come  to  zero. 
Figure  A-2  illustrates  the  effect  on  such  a  system  of  various  types  of 
viscosity.  The  first  profile  is  that  for  a  system  in  which  the  effective 
viscosity  has  been  removed  by  setting  f  =  -0.5,  in  the  second  case  no 
artificial  viscosity  has  been  applied,  while  in  the  third,  fourth  and 
fifth  profiles  an  f  type,  a  Richtmyer-Von  Neumann  type,  and  an  ac^  type 
viscosity  have  been  applied  respectively,  all  with  the  viscosity  coeffi¬ 
cient  set  equal  to  1 .  Notice  that  the  best  damping  is  obtained  with  the 
acQ  type  viscosity. 

The  effectiveness  of  the  artificial  viscosity  terms  depends  to  a 
large  extent  on  the  manner  in  which  they  are  expressed  in  the  difference 
equations.  This  is  true  not  only  in  the  PIC  method  but  also  in  other  dif¬ 
ferencing  schemes. 

Let  us  examine  the  effectiveness  of  Eqs.  (A-2J  when  they  are  modified 
to  include  artificial  viscosity.  The  obvious  way  to  achieve  this  modifica¬ 
tion,  without  destroying  the  conservativeness  of  the  system,  is  to  replace 
cell  pressures  in  those  equations  by  the  sum  of  the  equation  of  state  pres¬ 
sure  for  the  cell  plus  the  average  value  of  artificial  viscosity  at  the 
cell's  two  boundaries.  But,  by  reason  of  this  definition,  the  momentum 
equation  becomes 


(A-4) 
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tfiich 
id  re: 
:d  wi‘ 


The  defect  in  this  manner  of  differencing  is  apparent.  In  averaging 
the  dissipative  terms  over  such  a  wide  spatial  range,  the  sensitivity  of 
velocity  change  to  velocity  gradient  is  impaired.  This  results  in  too 
much  smoothing  in  the  neighborhood  of  shocks  and,  at  the  other  extreme, 
too  little  damping  when  perturbations  occur  in  an  otherwise  stagnant 
fluid. 

Nor  can  this  sensitivity  be  achieved  simply  by  writing  the  momentum 
equation  as 


(A-5) 


for  this  destroys  the  energy  conservation  of  the  equations. 

But  it  is  possible  to  achieve  difference  equations  which  are  conser¬ 
vative  and  yet  do  not  involve  interpolated  viscosity  terms.  To  do  this  it 
is  necessary  to  alter  the  derivation  of  the  Phase  I  equations.  Begin  with 
the  one  dimensional  momentum  and  energy  differential  equations  with  trans¬ 
port  terms  dropped  and  pressures  modified  to  include  artificial  viscosity. 
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But  now  rewrite  the  energy  equation  in  the  equivalent  form 
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Using  this  alternate  form  of  the  energy  equation,  we  write  the  equa¬ 
tions  in  the  difference  form 
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where  we  have  now  written  the  momentum  equation  in  the  form  of  Eq,  (A-5) 
With  the  difference  equations  thus  expressed,  energy  conservation  is  no 
longer  a  problem.  For  the  variation  in  total  energy  over  Phase  I  is 
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which  is  clearly  in  conservative  form. 

All  of  the  problems  described  in  this  report  made  use  of  the  comput¬ 
ing  scheme  described  by  Eqs.  (A-6). 
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